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ABSTRACT 

RX J1856. 5-3754 is the X-ray brightest among the nearby isolated neutron stars. Its X-ray spectrum is 
thermal, and is reproduced remarkably well by a black-body, but its interpretation has remained puzzling. 
One reason is that the source did not exhibit pulsations, and hence a magnetic field strength — vital input to 
atmosphere models — could not be estimated. Recently, however, very weak pulsations were discovered. Here, 
we analyze these in detail, using all available data from the XMM-Newton and Chandra X-ray observatories. 
From frequency measurements, we set a 2a upper limit to the frequency derivative of \v\ < 1.3 x 10~ 14 Hzs _1 . 
Trying possible phase-connected timing solutions, we find that one solution is far more likely than the others, 
and we infer a most probable value of v = (-5.98 ±0.14) x 10~ 16 Hzs _I . The inferred magnetic field strength 
is 1.5 x 10 13 G, comparable to what was found for similar neutron stars. From models, the field seems too 
strong to be consistent with the absence of spectral features for non-condensed atmospheres. It is sufficiently 
strong, however, that the surface could be condensed, but only if it is consists of heavy elements like iron. Our 
measurements imply a characteristic age of ~4Myr. This is longer than the cooling and kinematic ages, as was 
found for similar objects, but at almost a factor ten, the discrepancy is more extreme. A puzzle raised by our 
measurement is that the implied rotational energy loss rate of ^3 x 10 30 ergs _1 is orders of magnitude smaller 
than what was inferred from the Ha nebula surrounding the source. 

Subject headings: stars: individual: (RX J 1856.5-3754) — stars: neutron — X-rays: stars 



1. INTRODUCTION 

The nearby neutron star RX J 1856.5-3754 (hereafter 
J 1856) is the closest and brightest of the seven radio-quiet, 
isolated ne utron stars (I NS) discovered by ROSAT (for re - 
views, see lHaberlll2007l and Ivan Kerkwiik & Kaplanll2007h . 
The INSs have attracted much attention, in part because of 
the hope that the equation of state in their ultradense interiors 
could be constrained using their thermal emission. Progress 
has been slow, however, as adequate fits of the spectra appear 
to require somewhat contrived models, with the most success- 
ful one appealing to hydrogen layers of fin ely tuned thickness 
superposed on a condensed surface (e.g.. iMotch et al.1 120031: 
IHo et alJl2007h . A major hindrance is our ignorance of the 
surface magnetic field: without an observational constraint, 
models can consider too wide a range to be useful. 

We have been able to make some progress, using X- 
ray observations to determine coherent timing solutions and 
henc e estimates of the dipolar magn etic fields f or two of the 
INS (Kaplan & van Kerkwiill2005allR hereafter[KvK05lb). 
However, J 1856 has resisted such attempts, as its pulsations 
were so weak that the y were discovered only rec ently, in a 
long XMM obse rvation (iTiengo & Mereg hetti 2007, hereafter 
TM07). TM07 were unable to determine the spin-down rate 
of the source, finding only a limit of |z/| < 4 x 10~ 14 Hzs _I (at 
90% confidence). Here, from a detailed analysis of a larger 
amount of data, we infer a stronger constraint, and derive a 
most likely phase-connected solution. 

2. OBSERVATIONS 

We retrieved all observations of J 1 856 taken with the XMM- 
Newton (XMM) and Chandra X-ray Observatories (see Ta- 
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bleHJ, and reprocessed the data from XMM's European Pho- 
ton Imaging Cameras with PN and MOS detectors using the 
epchain and emchain pipelines in SAS version 7.1, and 
those from Chandra's High Resolution Cameras following 
standard threads in CIAO version 4.0. For the XMM data 
from 2007 Mar 14, the pipelines gave warnings about odd 
time jumps in both MOS and PN. We traced these to small 
sets of duplicated events in the observation data files (ODF), 
which we removed (our results do not depend on this, or on 
whether or not we include this observation). 

Given that the pulsations are so weak, we tried to optimize 
the number of source counts. For the PN imaging observa- 
tions, we decreased the default 150eV low-energy cutoff to 
100 eV, which leads to a 50% increase in source events. For 
both MOS and PN, we selected source events from a circular 
region of 37. "5 radius and with energies below 1 keV (where 
background flares are negligible). Following normal practice, 
we included only one and two-pixel events (patterns to 4) 
with no warning flags for PN, and single to triple events (pat- 
terns to 12) with the default flag mask for MOS (removing 
< 0.1% of the source counts). 

There are also two XMM timing observations. For the one 
with MOS1, we used default settings, and selected single to 
triple events from columns 315 to 329 with energies less than 
1 keV and the default flag mask. The timing observation with 
PN suffers from frequent bursts of low-energy noise events, 
and from some high-energy ones due to flares. We identified 
the former using epre ject, after which the source domi- 
nates the background down to 215 eV for single-pixel events 
and down to 430 eV for doubles, and up to 600 eV during 
flares and up to 800 eV otherwise. We selected source events 
using these criteria from columns 29 to 45. 

For the Chandra spectra, we e xtracted not only the zeroth 
order source events, as done by ITM071 but also the about 
two times more numerous diffracted events. For the former, 
we used a circle with radius of 3 "6, while for the latter we 
took rectangles in transmission grating angles (tg r ,tgd) de- 
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TABLE 1 

Log of Observations And Timing Measurements 
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(MJD) 
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2000 Mar 10. . . . 
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55 


36799 
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2001 Oct 8-15... 
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626 


320541 


2 


13 


0.9(2) 


52193.786169(4) 


0.1417396(3) 


2002 Apr 8 


Xmlti 


58 


664279 
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34 


1.01(17) 


52373.011797(2) 


0.1417403(17) 


2002 Aug 6 
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10 


9984 
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2002 Sep 3 
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55 


55632 


11 










2002 Sep 23 ... . 


Ifoc 
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14627 












2003 May 4 .... 
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50 


51544 
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2004 Apr 17 ... . 


Xpnti 


65 


377313 


9 


18 


1.0(2) 


53113.309266(3) 


0.141741(3) 


2004 Sep 24 ... . 


Xmmf 


71 


463180 




21 


1.0(2) 


53272.338862(3) 


0.1417401(19) 


2004 Nov 4 


I 


49 


47443 


14 










2005 Mar 23 ... . 


Xmllw 


35 


378678 




17 


0.9(2) 


53452.558481(3) 


0.141738(4) 


2005 Apr 15 ... . 


Xnopn 


9 


28357 












2005 Jun 10 . . . . 


I 


50 


65608 


2 










2005 Sep 24 ... . 


X 


35 


366372 




29 


1.3(2) 


53637.530095(3) 


0.141740(3) 


2005 Nov 13.... 


I 


48 


35586 


37 










2006 Mar 26 ... . 


X 


70 


772404 




16 


0.65(16) 


53821.055080(3) 


0.1417375(15) 


2006 Jun 1 


I 


46 


33270 


40 










2006 Oct 24 ... . 


X 


73 


794684 




52 


1.14(16) 


54032.442099(2) 


0.1417412(11) 


2006 Nov 18. ... 


I 


50 


31997 


40 










2007 Mar 14 ... . 


Xodf 


69 


734710 




55 


1.22(17) 


54174.246834(2) 


0.1417357(12) 


2007 Mar 25 ... . 


Xmlw 


41 


488457 




31 


1.1(2) 


54184.467408(3) 


0.141739(2) 


2007 Oct 4 


X 


70 


416891 




11 


0.7(2) 


54377.621398(4) 


0.1417405(19) 



NOTE. — Explanation of columns: (date) date of the observation; (ins) instrument code as given in note (a); (AT) 
time between first and last event; (N ev ) number of events extracted; (/j,) estimated fraction of events due to background 
(given only if > 1%); (Z^) power of the pulsations; (a) fractional amplitude; (TOA) time of maximum light closest to 
the mean of all event times; (u) frequency. Power and timing measurements were derived only for observations with 
/Vev > 2 X 10 5 . All uncertainties are lcr; those on TOA and v were derived from fits with a fixed to 0.96(1 —AW" ( see 
text). 

11 Instrument codes are: (S) Chandra High-Resolution Camera for Spectroscopy, used with the Low-Energy Transmission 
Grating; (I) Chandra High-Resolution Camera for Imaging, out of focus; (Ifoc) as (I), but in focus; (X) XMM European 
Photon Imaging Cameras, combining data from the PN and MOS detectors using small-window mode and thin filter; 
(Xmlti) as (X), but with MOS1 used in timing mode (for which = 8%); (Xpnti) PN used in timing mode, and both 
MOS used in full-frame mode and thus ignored; (Xmmf) for MOS, a portion was taken using the medium filters; (Xmllw) 
MOS1 used in large-window mode; (Xnopn) only the MOS were used; (Xodf) duplicate frames in the observation data 
files were removed; (Xmlw) both MOS used in large-window mode. 



fined by |tg d | < 0.°000531 (the default for extracting spectra), 
and 0.°12 < |tg r | < 0.°42 (corresponding to 21 < A < 73 A), 
where the source clearly dominates the background. 

Most of the Chandra imaging observations were not taken 
in focus, and we used ellipses to define regions where the 
source clearly dominated the background. We extracted 
events with pulse intensity in the range 1-220 only, since few 
source events had higher values. 

3. INCOHERENT ANALYSIS 

Given the fractional amplitude of only a ~ 1 % of the pulsa- 
tions in J1856, a 3-cr detection requires ~2(3/a) 2 ~ 2 x 10 
counts (for a single trial at a known frequency). Combining 
events from the PN and MOS cameras, ten out of eleven XMM 
observations have sufficient counts, while only the long spec- 
trum suffices among the Chandra obs ervations. For all of 
these, we computed Z\ power spectra (Buc cheri et al.|[l983h 
for the barycentered ev ent tim es in a narrow interval around 
the frequency found by TM07. As can be seen in Fig. Q] the 
pulsations are detected in all eleven observations. We then 
fitt ed the event t imes for each observation with a sinusoid us- 
ing |Caih| (|l979) minimization. The resulting fractional am- 
plitudes a, frequencies v, and inferred arrival times TOA are 
listed in Table |T| 

One sees that the different observations give similar am- 
plitudes. Correcting them for background and statistical bias 
(the expectation value is (a m ) = ([ao(l - /b)] 2 + 4/A f ev ) 1 ' 2 , 
where a m and ao are the measured and unbiased amplitudes, 
N ev the number of counts, and /b the fraction due to back- 



ground), and taking a weighted average, the best estimate of 
the amplitude is 0.96 ±0.06% (x 2 = 12.2 for 10 degrees of 
freedom [DOF]). 

To obtain the best measurements of the phases and frequen- 
cies, we refitted each observation holding the amplitude fixed 
at a; = 0.96(1 -/b)%. This yielded the same best-fit values, 
but uncertainties that were slightly increased for observations 
for which a m was high, and decreased for those for which a m 
was low. Analytically, this is expected: for small amplitudes, 
a is not covariant with v or cf>, and the uncertainties scale with 
(flm/fli)- 1 / 2 . In TableQ] we list uncertainties from this second 
fit. 

Fitting the frequencies gives a best-fit frequency deriva- 
tive v = (-5±4) X 10- 15 Hzs-', with x 2 = 13.7 for 9 DOF 
(see Table [2]). Clearly, one cannot exclude a constant fre- 
quency (which has \ 2 = 15 for 10 DOF). The 2-er upper 
limit is \v\ < 1.3 x 10~ 14 Hzs _1 (a factor of three improved 
over ITM07T, mostly due to the precise Chandra frequency). 
To verify our result, we also added the Z\ power spectra for 
a range of frequency derivatives; this led to the same best-fit 
values, and showed no significant other peaks (see Fig. [I). 

4. COHERENT ANALYSIS 

Since J 1856 has been observed so often, we attempted a co- 
herent analysis, using two different methods. First, we tried 
finding (v, v) combinations that were consistent with both the 
frequencies and the arrival times listed in Table Q] For this 
purpose, we generalised the method of lKvK05al to allow us to 
explore many possible cycle counts between observations (no 
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FIG . 1 .—Power spectra for RX Jl 856.5-3754. (left) Results for the long in- 
dividual data sets, near the pulse frequency identified by TM07 (dotted line). 
The ordinate is scaled such that the mid-point corresponds to a fractional 
pulsation amplitude a = (2Zj /JVev) 1 ' 2 = 1%. (top right) Incoherently com- 
bined power spectra for a range of period derivatives. The contours are at 
AZ? = 2.3, 6.17, 11.8 (or 1, 2, and 3<r for two parameters of interest). The 
circles represent the three highest peaks from the coherent analysis, (bot- 
tom right) Power for coherent combinations of the long observations (open 
circles) and all observations (filled circles). 

pair is close enough to give a unique starting solution), and to 
fit not just arrival times, but also frequencies. In more detail, 
using the long Chandra observation as a reference, we first 
calculated the number of cycles to the first XMM observation 
for the best-fit frequency and for v = 0. Next, we estimated 
an uncertainty on the cycle count using a ±5er range, with 
a determined from the measurement uncertainty on v and an 
assumed a priori uncertainty of 2 x 10~ 14 Hzs -1 on v. For 
each possible cycle count, we calculated a new, much more 
precise estimate of the frequency, and used it to estimate pos- 
sible cycle counts for the next observation. We iterated this 
process, deriving v as part of the fit for later iterations, until 
the fit clearly became bad (x 2 > 85, where including all ob- 
servations one has 19 DOF: 1 1 arrival times plus 1 1 frequen- 
cies minus three fit parameters), or until all observations were 
included. In the end, the best fits did not depend on which 
observation was used as an initial reference. 

The best trial, shown in Fig. [2] resulted in the frequency 
and frequency derivative listed in Table It has \ 2 = 25.7 
for 19 DOF. The fit seems reasonably unique: the two next- 
best possibilities have \ 2 = 49 and 53, then eight are found 
between \ 2 = 60 and 70, and thirteen between 70 and 80. Also 
ignoring frequency information, the best trial is superior: it 
has Xtoa = 10-^ f° r 8 DOF, while the next best has Xtoa = 24. 

As a second method of finding a coherent solution, we cal- 
culated Z\{v,v) for all data, coherently combining fourier 
spectra for individual observati ons for a range of frequency 
derivatives (Ra nsom et al.1120021) . This has the advantage that 
we do not have to decide a priori which peaks in the individ- 
ual power spectra are the correct ones (relevant especially for 
the long Chandra observation; see Fig . [T]> . For the combined 
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FIG. 2. — Phase and frequency residuals for RX J1856. 5-3754. (top) Fre- 
quency residuals relative to the best-fit quadratic model. ( center) Phase resid- 
uals relative to a linear (u = 0) model, (bottom) Phase residuals relative to the 
best-fit quadratic model. The filled circles represent XMM observations, and 
the open one the long Chandra observation. 



TABLE 2 

Timing Parameters for RX J 1856.5-3754 





Incoherent 


Coherent 


Quantity 


Value 


Value 


Dates (MJD) 


52,194-54,377 


?o(MJD) 


53000.0 


53000.000009(3) 


v (Hz) 


0.1417393(6) 


0.1417393685(5) 


i>(l(T 16 Hz s~')... 


-50(80) 


-5.98(14) 


TOArms (s) 




0.24 


\7l)OF 


13.7/9 


25.7/19 


P(s) 


7.05521(3) 


7.05520288(2) 


P(10" 14 ss _1 ).... 


20(40) 


2.97(7) 


£(10 30 ergs-') .. 


<70 


3.3 


Bdip (10 1 G) 


<7 


1.5 


Tchar (Myr) 


>0.17 


3.8 



NOTE. — T c har = P/2P is the characteristic age, as- 
suming an initial spin period Pq <^ P and a constant 



magnetic field; Bjip = 3.2 X 10 V PP is the magnetic 
field inferred assuming spin-down by dipole radiation; 
E = 10 45 /4547r 2 i/i> is the spin-down luminosity (with 
/ = 10 45 /45 g cm 2 the moment of inertia). Uncertainties 
quoted are twice the formal 1-cr uncertainties. For the in- 
coherent analysis, we used the 2<r lower limit of \u\ < 
1.3 X 10" l4 Hzs _1 to derive E, B, and r char . 

data, the highest peak has Z\ = 263 and occurs at the same 
(y,v) found above. This corresponds to a fractional ampli- 
tude a = (2Zf/7V tot ) 1 / 2 /(l - / b ) = 0.94%, consistent with what 
was found from the individual observations (here, the total 
number of events 7V to t = 6188356 and the background fraction 
/b = 1 .7%). The second and third-highest peaks have Z\ = 234 
and 230, and correspond to the second and ninth-best solution 
found using the trial-and-error method, with x 2 = 49 and 66. 
The changes in ordering arise because here we included the 
short observations, which causes the power in the main peak 
to increase by AZ 2 = 5, but that in most other peaks to de- 
crease (see Fig. [TJ. Folding the shorter observations on the 
different solutions confirms this conclusion. 

In order to verify the uniqueness of our solution, we ran 
1000 simulations in which we assumed our best solution and 
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created photon time series corresponding to each of the ob- 
servations (assuming a = 0.96(1 —fi,)%). We analyzed these 
in exactly the same way as the real observations. Among the 
simulations, in 983 out of 1000 cases the correct solution was 
recovered by the trial-and-error method on the long observa- 
tions, and in all 1000 using the Z\ power spectra on all data. 
Inspection of the mis-identification shows that, as expected, 
it is the addition of the information from the shorter observa- 
tions that causes the Z\ method to do better. 

The best-fit reduced \ 2 slightly exceeds unity, with 
X 2 /DOF = 1.4 (the second best solution has x 2 /DOF = 2.6). 
This could indicate a fundamental problem, but perhaps 
more likely reflects that for low-significance detections, out- 
liers in phase and frequency happen more often than ex- 
pected based on a normal distribution. Indeed, among our 
1000 simulations, we find that 211 have best solutions with 
X 2 > 25.7, somewhat more than the 140 expected for nor- 
mal distributions. Alternatively, some unmodeled phase var i- 
ations may be present, su ch as seen in other INS (KvK05a,b; 
Ivan Kerkwiik et alj|2007l) . 

5. RAMIFICATIONS 

Assuming the star is spinning down by magnetic dipole ra- 
diation, one can use the spin-down rate to infer a magnetic 
field strength, characteristic age, and spin-down luminosity 
(see Table [2]). We discuss the ramifications below, assum- 
ing the coherent solution is the correct one. We compare 
the results with those obtained for RX J0720.4-3125 and RX 
J 1308.6+2127 (J0720 and J 1308 hereafter). 

The inferred value of the magnetic field strength of 1.5 x 
10 13 G is similar to, but somewhat lower than the values 
of 2.4 and 3.4 x 10 13 G found for J0720 and J1308. This 
lower value might be consistent with the idea that the X- 
ray absorption features found in the other sources — but not 
in J 1856 — are due to proton cyclotron lines or transitions in 
neutral hydrogen, and that in the lower field of J 1856 these 
are shifted out of the observed band. The infe rred value, how- 
ever, is still somewhat high: calculations bv lHo et ail d2007l) 
suggest that for the approximate temperature of J 1856 and 
B = 1.5 x 10 13 G (and for a gravitational redshift zgr — 0.3), 
strong features due to bound-bound transitions should appear 
at energies of ~ 130eV (quantum number m = — » 1) and 
230 eV (m = — ► 2), but none are observed. One could ap- 
peal to elements other than hydrogen, but these generally have 
more bound transitions, thus exace rbating the situation (e.g., 
iPons et alJl200llMori & Holl2007l) . 

On the other hand, the inferred magnetic field may be con- 
sistent with the idea that the surface is condensed, and that 
this causes the black-body like spectrum . This depends on th e 
composition: from the recent work by Medin & Lail (120071) . 
it appears that for relatively light elements (based on calcu- 



lations for carbon and helium), J 1856 is far too hot for con- 
densation to be possible. For iron, however, the condensation 
temperature is relatively high: fcr con( j ~ 70 eV for a magnetic 
field of 1 .5 x 10 13 G. While we do not know the exact surface 
temperature of J 1856 bec ause of uncertai nties in redshift and 
color correction, the fits o flHo et all (l2007h have kToo ~ 40 eV. 
This corresponds to kT ~ 55 eV at the surface, suggesting that 
a condensed iron surface is possible. Indeed, this might also 
be the reason that in ea rly observations, J0 720 had a feature- 
less spectrum as well dPaerels et alj|200lb . For its field of 
2.4 x 10 13 G, iron could condense below kT con d ~ 1 lOeV, and 
its observed temperature kT^ ~ 85 eV corresponds to a sur- 
face temperature of ~ 1 lOeV. 

The characteristic age of 4 Myr we derive for J 1 856 is much 
larger than the kinematic age of 0.4 Myr inferre d assuming an 
origin in the Upper Scorpius OB association ( Walteri 120011 : 
Ivan Kerkwi ik & Kaplan 2007), and also greatly ex ceeds sim- 
ple es timates of the cooling age (~ 0.5 Myr; e.g., iPage et all 
2006). Longer characteristic ages — although by a factor of 
three rather than ten — were also found for J0720 and J 1308, 
which strengthens the suggestio n that this is a property com- 
mon to all isolated neutron stars dKvK05bt see lKvK05al for a 
discussion of possible causes). 

Like for J0720 and J 1308, the implied spin-down lumi- 
nosity £~ 3 x 10 30 ergs _1 is much smaller than the X-ray 
luminosity Lx — 3 x 10 32 erg s~ 1 (for a distance of 160pc; 
Ivan Kerkwi ik & Kaplan 2007), consistent with the lack of 
non-thermal emission. It is also, however, orders of mag- 
nitude lower than the independent estimate of E > 1 x 
10 33 (t//160pc) 3 ergs" 1 made by assuming that the Ha neb- 
ula associated with J 1856 is due to a bow shock, where the 
pressure from the pulsar wind matches the ram pressure from 
the interstellar medium dvan Kerkwiik & KulkarnilfeOOll) . In- 
deed, this discrepancy remains even if one considers just the 
incoherent anal ysis. The alternate model for th e Ha nebula 
considered by Ivan Kerkwiik & Kulkarnil d2001l) — that it was 
a moving i onization nebula (Blaes et alJI 1995D — was already 
rejected by iKaplan et alJ d2002l) . because the opening angle 
of the nebula's tail did not match observations for distances 
greater than lOOpc. Thus, our new measurement leaves the 
nature of the nebula an enigma. 
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